Level Curvature Distribution and the Structure of Eigenfunctions 

in Disordered Systems. 



C. Basu\ C.M.Canali^B, V. E. Kravtsovl'^ I. V. Yurkevich^ 
^International Center for Theoretical Physics, 
P.O.Box 586, 34100 Trieste, Italy. 
^ Dept. of Applied Physics, Chalmers University of Technology 
and Goteborg University, S-412 96 Goteborg, Sweden. 
^Landau Institute for Theoretical Physics, 
Kosygina str. 2, 117940 Moscow, Russia. 
^The University of Birmingham, School of Physics and Astronomy, 
Edgbaston, Birmingham B15 2TT, UK. 
(today) 



Abstract 

The level curvature distribution function is studied both analytically and nu- 
merically for the case of T-breaking perturbations over the orthogonal ensem- 
ble. The leading correction to the shape of the curvature distribution beyond 
the random matrix theory is calculated using the nonlinear supersymmetric 
sigma-model and compared to numerical simulations on the Anderson model. 

It is predicted analytically and confirmed numerically that the sign of 
the correction is different for T-breaking perturbations caused by a constant 
vector-potential equivalent to a phase twist in the boundary conditions, and 
those caused by a random magnetic field. 

In the former case it is shown using a nonperturbative approach that 
quasi-localized states in weakly disordered systems can cause the curvature 
distribution to be nonanalytic. In 2d systems the distribution function P{K) 
has a branching point aX K = Q that is related to the multifractality of the 
wave functions and thus should be a generic feature of all critical eigenstates. 
A relationship between the branching power and the multifractality exponent 
^2 is suggested. Evidence of the branch cut singularity is found in numerical 
simulations in 2d systems and at the Anderson transition point in 'id systems. 

I. INTRODUCTION. 

As first suggested by Edwards and ThoulesJil, the sensitivity of the spectrum {En} 
of disordered conductors to a small twist of phase (f) in the boundary conditions '^{x = 

is a powerful tool to probe the space structure of eigenfunctions 
and distinguish between the extended and the localized states. More precisely, the quantity 
Kn that is now referred to as the "level curvature" , was introduced in Ref. [ |l| in order to 
describe this sensitivity quantitatively: 
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(1-1) 
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where the mean level spacing A = {vL'^)~^ is related to the mean density of states v = {v{E)) 
and the size of the d-dimensional sample L. 

In complex or disordered quantum systems, the quantity Kn fluctuates over the ensemble 
of energy levels {E^} or, for a given level, over the ensemble of realizations of disorder. The 
typical width of the distribution of level curvatures P{K) is of the order of the dimension- 
less conductance g = D/{L'^A), where D is the diffusion coefficient. Thus studying this 
distribution one can study the Anderson transition from metal to insulator that takes place 
upon increasing the disorder. 

Recent work! has shown that the distribution of the fluctuating quantity is a partic- 
ular example of parametric level statistics, i.e. statistics of spectral responses of the system 
to a perturbation proportional to some parameter 0. 

A remarkable property that the parametric level statisticJl share with the usual level 
statisticsi'S, is that in a certain limit they are universal for all classically chaotic and dis- 
ordered systems and can be described by the random matrix theory (RMT) of Wigner 
and Dysoni'i. For disordered systems considered here this limit coincide Jl with g ^ oo. 
For chaotic systems the same role is played! by the ratio g = 7i/A, where 71 is the first 
nonzero mode in the spectrum of the Perron- Frobenius operator that describes the chaotic 
behavior of the corresponding classical system. In particular, for a time-reversal-invariant 
system without spin-dependent interactions (orthogonal ensemble) the distribution of level 
curvatures, Eq.( |l.l| ), was found&i in this limit to have the form: 

where {\K\)fi.mt = '^g is the average modulus of the level curvature. Further studyi has 



shown that the the form of the curvature distribution is still given by Eq.( |1.2| ) even when 
weak localization is taken into account; only the dimensionless conductance in the expression 
for {\K\) is changed appropriately. 

The form of Eq. ( |1.2|) is universal. It does not depend, e.g. on the details of the system and 
the perturbation. Its validity only requires the system to be T-invariant with = 1, and 
the perturbation to break this invariance. The underlying physics behind this universality is 
the basis invariance of RMT which is equivalent to the eigenf unctions of the physical system 
being structureless. 

Anderson localization apparently breaks the basis invariance. Thus at a sufficiently 
small value of g the universality of the spectral statistics should break down as well. In the 
strong localization limit one would expectS a logarithmically normal decay rather than the 
power-law tails of the Edwards-Thouless curvature distribution function. This is because 
the fluctuations of level curvature K oc e~^/^ can be viewed in this case as the consequence 
of the Gaussian fluctuations of the localization radius ^ for the exponentially localized wave 
functions. This picture is qualitatively confirmed by recent analytical0 and numericalElll^l 
calculations. 

In the present paper we present a comprehensive review of our recent analytical and 
numerical results on the correction to the level curvature distribution 6P{k) = P{k) — 
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Pwoik) when one approaches the Anderson transition point o = g*^ from the metal side 
g ^ 1. Some of the results discussed below are published 

It turns out that there are two completely different contributions to the correction 6P{k). 
One of them SPj-eg{k) is regular in the small parameter g~^ and can be obtained by a 
perturbative treatment^l of the nonzero spatial modes of the nonlinear supersymmetric 
sigma-modeli. The main result of this treatment0 is that the sign of the correction SPreg{k) 
depends on the topological nature of perturbation. It is different for the "global" Edwards- 
Thouless curvature, where the perturbation is represented by a global twist of the phase in 
the boundary conditions,lii and for the case in which the curvature is probed by a "local" 
T-breaking perturbation such as magnetic impurities or random magnetic fluxes. Below we 
present a numerical evidence of this fact. 

On top of the regular correction, there i^ also a nonperturbative in g~^ correction 
6Ps{k) which is proportional to exp{—l/g~^). The latter correction is due to the so called 
pre-localized statestOiHll, i.e. eigenstates with anomalously high peak(s) in the probability 
density |^'(r)p. 

There are reasons to consider the highly irregular, multifractal critical eigenstatesil as 
the result of a proliferation of such pre-localized states. This point of view is partly sup- 
ported by the observation that weakly localized states in the critical dimensionality d = 2 
also exhibit a (weakl multifractality, as can be shown by means of the same methods (renor- 
malization groupilo or space inhomogeneous saddle-point approximationEl) that were used 
to discover the pre-localized states responsible for the slow current relaxation in disordered 
conductor^HlH. 

This idea enables us to extend the results for 6Ps{k) obtained by the novel saddle-point 
approximationS'lli for 2d metals to the critical state at the Anderson transition in 2 -|- e 
dimensions. Thus we can explain the branching nonanalyticity at k = found numerically 
in0 for the 3d critical level curvature distribution function Pc{k). Furthermore, we suggest a 
relationship between the branching power and the exponent d2 describing generic multifractal 
critical statesil. This relationship fits well the numerical results and provides a link between 
the spectral statistics and statistics of wavefunctions near the Anderson transition. 

The paper is organized as follows. In Section II we describe the perturbative in g~^ 
approach for calculating 6Pj-eg{k). In Section III we generalize the instanton approximation 
of Ref. [ ^ for the problem of level curvature distribution and calculate the nonperturbative 
contribution 6Ps{k) for the metallic (weakly-localized) states in quasi-ld and 2d systems. 
In Section IV we extend the results of Section III to the critical states in 2 + e dimensions 
and derive the relationship between the branching nonanalyticity in Pc{k) and the fractal 
dimensionality In Section V the results of numerical simulation on the Anderson model 
are presented. Some open questions are discussed in the Conclusions. 



II. PERTURBATIVE CORRECTIONS TO P{K). 
A. Main results. 

A general approach to calculate the 1 /o-corrections using the nonlinear supersymmet- 
ric sigma-modelB has been suggested inli3 and applied to the distributions of different 
quantitiesi@0. It is based on a perturbative analysis of the nonzero diffusion modes 
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which are integrated out to produce 1/g corrections to the zero- mode supersymmetric sigma- 
modeli. The latter must then be handled exactly. 

Before going into the details of the calculations we would like to formulate the main 
results for dP^-egik) for the case of T-breaking perturbations over the orthogonal ensemble: 

2 — llA;^ + 2k'^ K 

where 



" {Trgy 5 {cfr ^ i -1, ' case ll " ^^'^^ 

Here q = {gi, ...qd}, where qi = 27ini, {rii = 0, ±1, ±2...) in the case of the periodic boundary 
conditions (for an unperturbed system) considered in this paper. 

A remarkable feature of Eqs.( p.l|) ,( |2.2|) ioi d < A is that the sign of the correction is 
different for global (case I) and local (case II) T-breaking perturbations. See Fig. |I[ 

The positive sign at small k in case I reflects the tendency towards a weaker spectral 
response to a change in the boundary conditions with decreasing g. It is related to long- 
range correlations in the wavefunctions that result in mesoscopic fluctuations of the matrix 
element of perturbation. The same long-range correlations cause mesoscopic fluctuations of 
the diffusion coefficient. 

The negative sign of the correction SPreg{k) for local perturbations [case II] is entirely 
due to the effect of the energy level statistics that lead to mesoscopic fluctuations of the 
density of states. The mesoscopic fluctuations of the matrix elements are suppressed in this 
case, because the effect of long-range correlations of the unperturbed wavefunctions is cut 
by the local nature of the perturbation. 

In order to illustrate the effect of the energy level statistics on the level curvature we 
invoke the expression for Kn in terms of the matrix element of the perturbation and 
the exact eigenvalues in the absence of perturbation En- 

Kn = 2Y. (2.3) 

From this expression we can clearly identify the two sources of the fluctuations in the level 
curvature. If the fluctuations of the matrix elements are suppressed, |KimP can be replaced 
by a constant. Then there is only one source of fluctuations left, that is the energy level 
statistics. Notice that implicitly this contains the effect of the statistics of the eigenf unctions 
as well, and in particular, their long-range correlations. Upon decreasing g the system of 
energy levels becomes less and less correlated. In the extreme limit of uncorrelated levels the 
distribution of curvatures is knownS to become of the Cauchy-Lorentz form Pciik) = l/7r(H- 
k"^). Since Pcl(O) = 1/vr and PwoiS^) = 1/2 we conclude that the effect of softening the 
energy level correlations on the shape of the curvature distribution is such that 5P(0) < 0, 
in full agreement with the results for the case II. 

The principal result of this Section is that Thouless relationship of proportionality 
{\K\) oc g breaks down beyond RMT. The ratio r{g) = {\K\) /2g increases above its RMT 
value r = 1, the correction being equal to: 
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2g (.,)2l.jqT I I, case II ' ^''"^ 



B. Functional representation for P{K). 



We know describe how Eqs. ( ^yiD , (p.2|) have been obtained. Let us consider a disordered 



mesoscopic d-dimensional system with a random white-noise impurity potential V^(r) per- 
turbed by a small vector-potential (f) /L. It is described by the microscopic Hamiltonian of 
the form: 

In case I we assume the sample to be closed into a ring geometry pierced by a small static 
magnetic flux $. Then the effect of the perturbation is equivalent to a twist of the boundary 
conditions generating the phase = 27r$/$o ($o = hc/e is the flux quantum). As a vector 
potential we will take a constant vector of the form /L = {(f)/L) n = (0/L){l, 0, ..0}. 

In the "local" case II we consider (r) to be a random 5-correlated vector-potential: 

(0,(r) 0^(r')) = v^(f)'' 5{r - r') = (2.6) 

The parameter is introduced in this way in order to keep {\K\)jiMT = '^9 the same as to 
the "global" curvature (case I). 

The curvature distribution function 

P{K)=^l^5{K- K,,) 6 {E - E^)^ (2.7) 

can be expressed in terms the two-level parametric correlation function R{uj,(j)) 

R (uj, 0) = (z/ {E + uj,<p)u {E, = 0)) (2.8) 



in a form similar to one derived in Ref. [ ^ for the distribution of level velocities: 

.2 . 1 



P (K) = Ihn [uj = -irA0^ 0j . (2.9) 

Indeed, using the exact expression for the fluctuating density of states ^'{E, 0) = 
L'"^ En ^{E - En{(p)) we have: 

R{uj, 0) = Y,{6{uj + E„(0) - E„(0)) 6{E - E^{0))). (2.10) 

m,n 

Because of the level repulsion, in the limit u ^ and 0^0 only terms with n = m 
contribute to the sum Eq.( |2.10D . On the other hand, with a perturbation which is odd 



5 



under time reversal, the T-invariant level energy ^ EniO) + \K„A(f must be eve ni 

in 0. Then choosing uj = ^KA 0^ we immediately arrive at Eq.( p.9|) . 

The two- level correlation function R{uj, 0) can be represented in the form of a functional 
integral using the Efetov's supersymmetry approach. A straightforward application of the 
results of Ref. [ || and Eq.( |2.9| ) leads to: 



64 dJidJ2 



Jl = J2=0 

where for the case I the functional F[Q] takes the form: 

, 2 



Z = \iin!^(f)'^^J DQexp{~F[Q])Y (2.11) 



F[Q] = -^JyStr\{LVQ + z <I>,Q 



2ik(j)\AQ) 



yStriJQ), 



(2.12) 



A similar representation for P{K) has been used in Ref. 
In Eq.( |2.12D we have introduced the notationS 



k=K/(2g), 



0=0 P+r, J = k [JiP+ + , 

and 

n = ^(l + A), P- = ^(l-A). 

The coordinate dependent 8x8 supermatrices Q (r) are parametrized as Q = T~^AT, 
where T belongs to a graded coset space UOSP (2, 2|4) / 
UOSP (2|2)®f/05P(2|2)i. 

Other matrices are specified as follows: 

A = diag {h, h, -h, -h)n^A ' ^ = '^^^5' (ts, ^3, 0, 0)^_^ , 

k = diag{l2,-hJ2,-h)R^A^ ^"3 = c^^^ctfi' (1, -1) , h = diag {1,1) . 

Above we imply the following hierarchy of blocks of supermatrices: retarded- advanced 
{R — A) blocks, boson-fermion [B — F) blocks, and blocks corresponding to time rever- 
sal. 

In the case II the linear in term in Eq.( p.l2| ) is absent but otherwise the functional F[Q] 
is the same provided that is introduced as in Eq. ( |2.6|) . A similar functional F[Q] appearsQ 
if one considers a small concentration of magnetic impurities as a perturbation. In both 
cases the structure of the "covariant derivative" VQ = VQ + 0, Q , which implies a sort 
of global gauge invariance, is broken down. 

It is important in deriving the functional F[Q] for the case II that the correlation radius 
of the random vector-potential is much smaller than the elastic scattering length. In this case 
the averaging over (r) should be done before switching to Q- variables that are assumed 
to be slowly varying in space. In the opposite limit of large correlation radius, one should 
average e"'^'*^' over (r) and arrive at a much more complicated functional. 
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C. Perturbative treatment of nonzero modes. 



The representation given in Eqs.( p.li| )-( p.l2| ), in terms of the field Q{r), contains all the 
spatial diffusion modes 7^ = (D/L^)q^. However, in doing the limit — in Eq.( p.llj ) the 
main role is played by the zero mode that corresponds to q = 0. At = this mode is 
gapless and thus it does not cost any energy no matter how large are the components of 
the field Q in the noncompact boson-boson sector!. It is the arbitrarily large amplitudes 
of the zero mode components of the field Q that compensate the infinitesimal parameter </> 
in Eqs.( |2.11] )-( |2.12| ) and lead to a finite result for P{K). Thus the space independent zero 



mode Qo must be considered nonperturbatively. 

In the limit (7 — 00 all the nonzero modes can be neglected!, and one arrived at the 



RMT result, Eq.(L2). For finite 1/g the nonzero modes should be also taken into account. 
However, all the nonzero modes can be treated perturbatively for (7 3> 1 leading to some 
corrections to the zero-mode action. In order to obtain these corrections we have to separate 
the zero modes from all other modes and then integrate over all the nonzero modes using a 
certain perturbative scheme. 

Following the method suggested by Kravtsov and MirlinS we decompose matrices Q (r) 
as follows: 

l + W/2 

Q{r) = T,-'Q{v)To, g = A-X^ (2.13) 

1 - W/2 

where Tq describes the zero mode and W{r) = I]q^o Wq^e^*^^ does not contain the zero mode 
at all. 

As has been already noticed, the main contribution to the functional integral comes from 
the zero mode. The zero-mode approximation, Q = Qo = Tq^^Tq is known to be equivalent 
to RMT§I. To go beyond RMT we integrate perturbatively over W{r) to obtain the effective 
zero-mode action F'^-l'-^ [Qq] as follows: 

F'" [Qo] = -\nl DQ-J[q] exp {-F [Qo, Q] } , 



where F Qo,Q is obtained from F[Q] by substituting the decomposition of Eq.( p.l3|) and 

J[Q] is the Jacobian of the corresponding nonlinear transformation. 

This scheme is implemented in Appendix ^. As a result we have the effective action 
expanded up to the second order in 1/g: 

peff ^ peff ^ peff ^ peff ^ ^e//_ ^2. 14) 



The first term in Eq. ( |2.14| ) Fq^^ [Qo\ is nothing but the zero-mode action responsible for the 
i?MT-resuhB: 

F^ff = \a-ikl3-Y.Jpap, (2.15) 

where p = ± labels the {R — A) blocks and 

a = STr(j)Qo)\ p = STr (^^Qo) , = STr[kQl'']. (2.16) 
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The next term Fi^^^ [Qq] is the first order (weak locahzation) correction obtained by Fyodorov 
and Sommersi: 

Fff = -^n,.a, (2.17) 

where 

The first order correction F^^^ leads to a renormahzation of the coefficient in front of a 
in Fq^^ and can be absorbed in the dimensionless conductance g. It can be checked that 
the renormalized coefficient g is exactly the conductance with weak localization corrections 
taken into account. In what follows we shall assume g to he a renormalized conductance 
and we shall omit Fi^-^ . 

We now consider the higher-order term F2'^-^ [Qo]'- 



A d d{d + 2)J \2 d d{d + 2 

+ (^k^p^ + ik{l - A/d) a(3 + a), 

and a higher-order source- induced contribution jFj [Qo] : 

Fff = n4 ■ I Q(l - A/d) a-ikp^Y.Jp<'p-X\.Jp ^p} • (2-20) 

The terms in F'^^f containing the factors 1/d and l/d{d+ 2) originate, after the angular 
integration over q, from the gradient term linear in 0, Str[(f)'VQ Q], in Eq. (|2.12| ). Such 
term is present only in the "global" case I. For this reason in the "local" case II all the 
(i-dependent terms in F'^-^-l' should be omitted. 

Differentiating the partition function Z with respect to the sources we arrive at the 
expression for the level curvature distribution function P (k), where k = K/{2g): 



P (k) = hm -^3? [dQ,-^ {a, (3) ■ exp " 



-0 327rg 



-+ikp 



StrkQl^StrkQ 



' 



$ (a, /3) = 1 - F^^^ + n4 (1 + 2ikp - (1 - A/d) a) . (2.21) 

A remarkable feature of the perturbation theory is that the function $ (a, (3) is a polinomial 
in a and This property is due to the fact that long-trace vertices that appear after the 
perturbative integration over W{r) factorize into the product of short-trace vertices a,j3, 
and (Tp. 

Then the correction 6P{K) can be represented as a finite order differential operator acting 
on the RMT distribution function Pwd{K) that corresponds to $ (a,/?) = 1 in Eq. ( |2.21| ): 
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where 



0- ^Pwoik/a) 



(2.23) 



lim 



-0 327r^ 



^ / DQ exp 



a + tkn 

2 ^. 



StrkQl^StrkQl 



Using the identities (see Eq. [A16| in Appendix |^) that relate the derivatives of Pa{k) with 
respect to k to those with respect to a, we can rewrite 5P{k) in the following form: 



5P {k) = U4 



(9 16 36 \ d\ 
+ I 2 T + did + 2) jYa]^'' ^^^^^=' ■ 



(2.24) 



Note that by the definition given in Eq.( p.23| ), the function Pa{k) obeys two normalization 
conditions: 



JdkPa{k) = l; Jdk\k\Pa{k) 



a. 



(2.25) 



Using the first of these conditions we immediately conclude that the cancellation of the 
terms proportional to Pa{k) in Eq. (|2.24|) ensures the conservation of the normalization 
Jdk6P{K) = 0. 

Next we note that the terms proportional to the first derivative can be absorbed into the 
function Pa{k): Pa=i{k) + 6{d/da) Pa=i{k) ^ Pa=i+s{k). In doing so we observe with the 
help of the second normalization condition in Eq. (|2.25|) that = {\K\) /2g = 1 + 6. Thus 
the terms with the first derivative in Eq.( p.24D result in a shift in the average (|-ft'|): 



(2.26) 



By redefining the k = K/ where {\K\) is the average of absolute value of the level 

curvature, one can cancel the terms with the first derivative in Eq.( 2.24 ). All what is 
left is the term with the second derivative which describes the change in the shape of the 
distribution function. The final result of these tedious calculations is very simple: 



5P{k) 



1 n4 



da? 



{a ^PwD{k/a) 



(2.27) 



a=l 



This equation means that the RMT curvature distribution PwD^k) plays the role of the 
generating function for its own corrections. 



III. THE SIGNATURE OF THE PRE-LOCALIZED STATES IN THE LEVEL 

CURVATURE DISTRIBUTION. 

It has been known for quite a while that the relaxation of current and the local density 
of states (DOS) in disordered conductors exhibit an anomaly even in the weak-localization 
regime. Namely, it has been shown in Ref. [ that there exists a small (but not ex- 



ponentially small) probability of finding a current relaxation time or a local DOS that is 
much larger than the corresponding mean values. These anomalies have been attributed to 
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quasi-localized (or pre-localized) statesEjE2l, that is, states with an anomalously large peak 
in |\E'(r)p at some point r = Tq. 

Very recently the problem of current relaxation in disordered conductors has been 
reconsiderecfiij by an elegant instanton approximation^ applied to the supersymmetric 
version of the nonlinear sigma-modeli. In these papers the main result of the previous 
wor has been confirmed for 2d systems. However, the new method was able to describe 
some unknown regimes of current relaxation and to set correct limits of validity for the 
regimes found earlier. Later the same ideaS has been appliedEl to find directly the distri- 
bution of |\E'(r)p and the distribution of local densities of states^. Thus the existence of 
quasi-localized states has been proved and the corresponding configuration of the random 
impurity potential has been foundii'0. 

The main idea of Edwards and Thoules^l is that it is possible to distinguish between 
localized and extended states by analyzing the sensitivity of the spectrum to a twist of 
the phase in the boundary conditions. This sensitivity is significant only for states with a 
localization radius larger than the sample size L and negligible for strongly localized states. It 
is clear that the existence of the pre-localized states should lead to an enhancement in P{K) 
at small K. In low dimensional systems d = 1,2 where the pre-localized states correspond to 
localized states with an anomalously small localization radius, one may expect a singularity 
in P{K) at i^' = 0. In 3d metal the typical pre-localized state looks like a sharp peak in 
|\&(r)p on top of the extended background |\E'(r)p ~ L~'^ = const. The level curvature 
that corresponds to such a state does not vanish but only slightly decreases. Thus the pre- 
localized states in 3d in the weak localization regime should have much weaker effect on the 
level curvature distribution. 



A. Instanton approximation. 

In order to check these predictions we consider, instead of P{K) at small K, its Fourier- 
transform P(A) = JdKP{K) e~'^^ at A > 1. 

It is easy to see that for both the RMT result Eq. (|1.2| ) and the regular correction Eq.( p.l|) 
the function P{X) vanishes exponentially for A 3> 1. In what follows we will seek for slowly 
decreasing contributions. Support for the existence of such contributions can be gained 
by noticing that in the functional representation of P{K), Eq. ( p.llD , Eq.( |2.12| ), the level 
curvature k plays the same role as the frequency in the problem of the current relaxation in 
disordered conductorJ^l, thus A being analogous to time. Therefore, one may expect long 
nonexponential tails in P{\) in analogy with those present in the current relaxation function 
J(t). However, it is far from clear that two problems are equivalent, since the boundary 
conditions are different and the nonlinear sigma- model in Eq. (|2.12| ) contains additional terms 
that describe the T-breaking perturbation. 

The main idea of Ref . [ ^ , that we will exploit here, is that at large A the configurations 



of the field Q{r) that are space independent or slowly varying in space, are energetically 
unfavorable. In contrast, essentially space-dependent configurations in the vicinity of the 
classical (instanton) solution Qinsi^^) that minimizes the action F[Q] appear to be energeti- 
cally advantageous. At large g the fluctuations around this solution are expected to be small 
and one arrives at: 
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//•H-oo 
VQ rfirA[Qi„,;0]e-™"=l+'^^>. (3.1) 
J-oo 

where A[Qins;(p] is a pre-exponential factor including the effect of fluctuations around the 
instanton solution. 

The Grassmann variables in the action, Eq. (|2.12| ) can lead only to a renormalization of 
the pre-exponential factor A in Eg. ( p. ID , since the integration over these variables is equiva- 
lent to a differentiation. Thus with exponential accuracy we can neglect all the Grassmann 
variables in the Efetov's parametrization for Q{r). Next, a finite contribution to S{X) in 
the limit 0-^0 comes only from the infinitely large boson-boson components of the field 
Q{r). Therefore we consider only the leading terms in the noncompact angles 6i and 62 in 
the Efetov's parametrization for the orthogonal ensemble: 



Q = V-^HV, (3.2) 



where 



with Pb = {k + l)/2, 



cosh 9b = cosh 9i cosh 62 + o"x sinh 9i sinh 62, (3.4) 
sinh 9b = sinh 9i cosh 6*2 + cr^ cosh 9i sinh 6*2, 

and 

V = e'"""^ ®P+®Pb + e^^"^ ® P_ ® Pb- (3.5) 

The field Qinsij) must obey periodic boundary conditions, since the twist of phase is 
taken into account explicitly in the action This means that ^1,2 (r) and the functions 

exp[(y9(r)], exp[x(r)] should obey periodic boundary conditions. 
In this way we obtain F[Qins\ = L~'^ J f[Qins] d'^^, where: 

f[Q] = ^g[{d9^f + id9^f]+ (3.6) 



TC 

2' 



{d^ - 0n)2 + {dxf (cosh 9+ cosh ^_ - 1) - 



tlx 

— 2dx {dip — 0n) sinh 6*+ sinh 9_ — '^^^'^ (cosh 6*+ + cosh 9_ ) . 

Here 9± = {9^ ± 92)/2 with 9i,2 > and 9 = {d/dxa}. 

We will look for a minimum of the functional F[Qins] + iKX that corresponds to: 

9_ = dx = 0, 9+ = 9e [0, +00]. (3.7) 

By varying the functional F[Qins] + iKX over 9, ip and k we find: 

92^ + 0^[K-(9t;-n)2] sinh = 0, (3.8) 
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d[(dv -n){cosh9 -1)] = 0, (3.9) 

and 

j(f)'^ J{coshe + l)d'^p = X, (3.10) 
where k = iK/2g, d'^p = ^ and = 

Eqs.( |3.8| )-( |3.10| ) correspond to the global case I. As usual, in the local case II the terms 
linear in n = {1, 0, 0, ...} are absent. 

The limit ^ is done simply by absorbing 0^ into 6. We introduce 6 = + ln0^. 
Then in the limit ^ we have sinh6' ^ cosh6' = ie^-^ ^nd Eqs.(|3^),(F^, (|3lOD take 
the form: 

d^e + ^[K-{dv-nf]e' = 0, (3.11) 

d[{dv -n)e^] =0, (3.12) 

I [ e'd''p = X. (3.13) 
8 J 

where now 9 G [—oo, +oo]. 

Using Eqs. (PTB|) , (|3 . 7|) , ( PTTTD and the periodicity of the function ^^(r) we find: 



S{\) = jgj {def d'^p + 2gK\. (3.14) 

We note that drops from the problem only if we assume a topologically trivial solution 
corresponding to periodic boundary conditions being imposed on v{v). Otherwise appears 
in the boundary condition for v = (f{r)/(j) which is periodic modulus 27r/0 and thus is ill 
defined in the limit — 0. In what follows we consider only such a topologically trivial 
solution. 

One can solve Eq. (|3.12|) : 



{dv -n) = [V X A]e-\ (3.15) 

where [V x A] = const in Id and is the curl of an arbitrary vector function A(r) in higher 
dimensions. Below we consider only the simplest solution that corresponds to [V x A] = 
—n/N = const. 

Let us consider first the local case II. Doing the space integration of Eq. (|3.15|) which in 
this case does not contain the term proportional to n, and using periodic boundary condi- 
tions for f (r) one immediately arrives at [V x A] = dv = 0. Then the same procedure with 
Eq.( p.ll| ) leads to the conclusion that the only solution for 9 that obeys periodic boundary 
conditions, is space- independent and exists only for k = = 1. The corresponding ac- 
tion is S{X) = 2gX. Thus the instanton approximation in the local case II gives only an 
exponentially small tail -P(A) oc e""^^^ that has been already obtained by the perturbative 
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approach. We conclude that for case II the analogy with the problem of current relaxation 
appears to be wrong. 

Now consider the global case I. Integrating Eq. ( |3.15| ) over space and using the periodicity 
of f (r) gives: 



N= I e-U'^p. (3.16) 



Substituting Eqs. (|3.15| ),( |3.16| ) into Eq. ( p.ll| ) we finally arrive at: 



d^e + ^ = d'e + - - e-' = 0. (3.17) 

It appears that the global nature of perturbation and the corresponding linear in n term 
in Eq. ( p.l5| ) leads to a term proportional to in Eq. ( p.l7| ) that builds a second "wall" in 
the effective "potential" U{6) and makes it possible for periodic solutions ("oscillations") to 
exist. 

Eq.( |3.17| ) takes a more symmetric form if we make a shift 6 = u — (, where: 

^oshC=(K+^) sinhC=(K-^) J^. (3.18) 



Finally we have the system of equations: 

d\ + -f^smhu = 0, (3.19) 



l^ = l' I e-^d'^p, (3.20) 



Je^'d'^p, (3.21) 



where 7^ = y/R/N. 

Solving these equations for a hyper-cubic sample —1/2 < < 1/2 with periodic bound- 
ary conditions one finds u{r, A), A^(A) and 7(A) which enter the instanton action S{X): 



SW = -9 I {duY d'^p + 2g^'N'\. (3.22) 



B. Non-exponential tails of P{\) in low-dimensional systems. 

We will see below that for large A the parameter 7 is small. For 7^1 the term 7^ sinh u 
is very small unless sinhu is exponentially large. This means that we can approximate 
7^sinhM ^ Y^'"' sign{u). Thus Eq. (|3.19|) is replaced by the Liouville equation : 

0,2 

d'^u + -^el"l sign{u) = 0. (3.23) 
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The generic solution to the Liouville equation in the low- dimensional case c? = 1, 2 is given 
in terms of two arbitrary functions f{w) and v{w) of the complex variable w = i'^zj\f% with 
z = X \ iy: 

eM = 2/'H^'(-"^) (3 24) 

where /' = df /dw and v = dv/dw. We need the r.h.s. of Eq. (p.24|) to be real positive. This 
can be done by the choice: 

r{w)v{-w*) = 1. (3.25) 
Then we have the solution in terms of one function only F{z) = f{i'yz): 

e'"' = \lL.. - (3-26) 



1. Quasi-ld Case 



Choosing F{z) = e^^^^ in Eq. (p.26|) , one has a quasi-lc? solution to the Liouville equation 



that depends only on one coordinate x: 

gM ^ (3.27) 

7^ cosh^(A;x + h) ' 

where k and h are real constants. 

The solution on a ring — | < a; < | is constructed by reflecting anti-symmetrically the 
positive solution with 6 = around the points x = ±| [see Fig. ^. The second constant k 
is found from the condition m(±1/4) = 0: 

4fc2 = 72 cosh2(fc/4), k ^ ln(l/7^). (3.28) 

The anti-symmetry of the solution immediately leads to the identity: 

e"cix = / e-^'dx. (3.29) 



Since the function u ~ ln(l/72) is large everywhere except in the vicinity of its zeros at 

\x\ 



i we have: 



8k 8k 



I ^ e"rfx = — tanh(fc/4) ^ — . (3.30) 

j_i Y 7 

Next we calculate the integral: 

J\\ufdx = 2 p~\ufdx = Ae - IQk ^ 4k^. (3.31) 



2 4 
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Then from Eqs.( p.20| ),( p.21|) we have: 

N = 



1 

8k' 



(3.32) 



7' = (3.33) 

Finally using Eq. (|3]2g) , (g;2|) and (|333|) we arrive at: 

5(A) = 7r^A;2 ^ TT^ln^A. (3.34) 
Thus the characteristic function in a quasi- Irf systems is: 



P(A) = Aexp 



-^In^A 
2 



91 = 27^9, (3.35) 



The above result holds within the domain of validity of the nonlinear sigma-model Eq. ( p.l2| ). 
This model and hence the saddle-point equations work only for sufficiently slow varying fields 
Q(r), namely \du\ < L/l, where / is the elastic scattering length. It follows immediately 
from Eqs. (|3.27|) , (|3.28|) that the above result is valid for 1 ^ A ^ exp{L/l). 

The logarithmically-normal tail in P{X) described by Eq.( |3.35| ) is exactly of the same 
functional form as the current relaxation function I[t) in Ref. [ ^ for the orthogonal en- 
semble. 



2. 2d Case 

In full analogy with the quasi-l(i case, we construct a double-periodic solution to the Li- 
ouville equation on a torus —^<x,y<^hy reflection. We consider a positive solution u{z) 
inside the square Q with vertices at z = ±1/2, ±i/2 and then continue it anti-symmetrically 
about a side of the square in any quarter of the sample \^z\ < 1/2, \Qz\ < 1/2. The defini- 
tion domain of the solution with its sign is drawn in Fig. ^. By construction, the symmetry 
relationship Eq. ( p.29|) is valid for such a 2c? solution too. 



The procedure of finding the solution is described in the Appendix y. We note that for 
our purposes we need only the solution for = r -C 1. It is rotationally invariant and has 



the form: 



_ 16&(fc - 1) V^^-^ 



where 

6 = 16( — 1 {k-lf. (3.37) 



' 9 \ fc 



Note that the solution Eq.( p.36| ) can be immediately obtained from the radial Liouville 



equation, with k and b being two constants of integration. The requirement of periodicity 
of u{z) helps to establish a connection, Eq. (p.37|) , between these constants. 
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The remaining constant k is found in a standard way from the requirement of convergence 
of J{du)'^dxdy in the action S{\): 



Thus we immediately find: 
and 



2k 



(du) 



-4b{k-r 
--2, b = 47r^ 



„2fc-3 



(3.38) 

(3.39) 
(3.40) 



(7^ + br'^y 

Because of the symmetry of u{z), the integral jQ{du)'^(Pp over the square Q is exactly 
one-half of the total integral over the period (over the sample) ^{duY (P p. It diverges 
logarithmically at r S> 7, and we arrive at: 

j {duf <fp = 327r In , (3.41) 

where C = can be found from an exact solution in the region l^l ~ 1. 

The result is almost independent of b at small 7 and is essentially determined by the 
logarithmic solution of the Poisson equation that follows from Eq.( p.l9D at 7 = 0. 

Now let us calculate the integrals in the self-consistency equations. For symmetry reasons 
we have: 

IGvr {k — 1) IGvr 



J e" dxdy = J e ^ dxdy 



(3-42) 

Most of the contribution to these integrals comes form the small r ~ 7 -C 1 and the result 
is independent of b. 

Now we are in the position to calculate the constants 7 and that enter the instanton 
action, Eg. (13^2^) . They are given by Eqs.(|]20D,(|3;2l|),(|j|): 



3^4 



N 



1 



8% 



167r' ' A 
Then the final expression for the instanton action in 2d reads: 



5(A) = A-K'g 



In 



- 1 



(3.43) 



(3.44) 



Accordingly, the characteristic function P{\) turns out to have a power-law asymptotic 
behavior at large A 3> 1: 

2S2 



P(A) = A 



92 = 27T g, 



(3.45) 



where c = 8e. 

Few notes should be made on the validity of the result Eq. (|3.45|) . Firstly, the above 
instanton approximation with the action S'(A) logarithmic in A is only justified when g 3> 
1, since the pre-exponential factor A could also be a power-law function of A but with 
an exponent of order 1. Secondly, the nonlinear sigma-model and hence the saddle-point 
equations work only for \du\ < L/l, where / is the elastic scattering length. It follows 
immediately from Eq. (|3.40|) , (|3.43| ) that the above result is valid for 1 ^ A ^ [L/l)^. 
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C. Nonanalyticity of the level curvature distribution. 

In this section we show that the slowly decreasing tails in the characteristic function P{X) 
at A 3> 1 given by Eqs.( |3.35| ), ( |3.45| ) result in a nonanalytic behavior of P(yK) at K = 0. As 
usual, true nonanalyticity arises only in the thermodynamic limit L/l ^ cxx3, since only in 
this limit the tails extend to infinity. For any finite L/l the function P{K) is still analytic 
at = but the region of the regular behavior of P{K) shrinks to zero with increasing 
L/l. Below the limit L/l — >■ cxo is assumed. 

Let us consider the quasi-lrf case first. In this case all derivatives of P{K) are finite at 
K = {]■ 

/ + 00 _ rjX 

^_iyy2n p(^) ^ ^^p[(2n + if /2g,]. (3.46) 
-oo ZTX 

Yet the function P{K) is nonanalytical at K=0, since the Taylor series P{K) = 

J2n ^(2n)''!°^ radius of convergence because of the very fast growth of P^'^'^\0) 
with n. 

The singularity at = is much stronger in 2d case. In this case all derivatives p(2n)(o) 
with 2n + 1 > 2g2 are proportional to (L/Z)^^^"'"''^"^^^) and diverge in the thermodynamic 
limit. Let us define m as an integer obeying the inequality of \g2 — ml < |. Then the the 
expansion of P{K) at small K has the form: 

P{K) = Co + CiK^ + ...c^_iir2(™-i) + c^iT^'"-"'" + oiK^""), (3.47) 
where < < 2 is given by: 

a„ = {2n + 1) - 2(72. (3.48) 

In the 3d metal case we failed to find a solution to the saddle-point problem that would 
lead to a finite action S{X) in the thermodynamic limit. This means that the characteristic 
function P(A) has only regular corrections at (7 ^ 1 and thus decays exponentially for A ^ 1. 

IV. NON-ANALYTICITY OF P{K) AT K = AND MULTIFRACTALITY OF 

EIGENFUNCTIONS. 

From the results of the previous section we see that the strength of the singularity of 
P{K) aX K = depends on the dimensionality in a nonmonotonic way. The singularity 
is very weak in a quasi-ld metal; it reaches a maximum in a 2d metal where P{K) has a 
branch cut; it disappears in a 3d metal where the level curvature distribution is analytic. 

Such a behavior is related to the fact that d = 2 is the low critical dimension for the 
Anderson transition, and the wavefunctions in the 2d weak-localization regime share some 
features of the critical wavefunctions at the Anderson transition in higher dimensions. 

A. P{K) at the Anderson transition in 2 + e dimensions. 

The usual way to describe the critical state near the Anderson transition is the {d—2) = e- 
expansion. To this end one considers the quantity of interest in a 2d system with g2 ^ 1 
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and then replaces q ? b y the critical conductance = 1/ (c? — 2) which is the fixed point of 
the scaling equationi 



3334 



— ^ = rf-2 -- + - . 4.1 
dlnL Qd \gij 

For the orthogonal ensemble in d = 2 + e dimensions we find to the leading order in e -C 1: 

P,(A)ocQ)\ (4.2) 

Note that at the critical point, the conductance is exactly size- independent, and one can 
consider the thermodynamic limit L — oo without tuning other parameters in order to keep 
9*d fixed. 

Thus one can define the critical exponent n that determines the power-law tail of the 
critical characteristic function: 

P,(A)ocA-^ /i = - + o(l). (4.3) 

e 

If we set e = 1 in the above equations we find /i ~ 2. Then it follows from Eq. (|3.47|) that 
already the second derivative of P{K) at if = is divergent: 

Pc(if) = co-ci|if|2-", a = 3-/i. (4.4) 



B. Exponent fi and multifractality. 

Unfortunately it is known that the accuracy of the d — 2 = e expansion is quite poor 
and insufficient for a precise determination of the critical exponents. In this situation one 
can try to find relationships between different critical exponents rather than try to evaluate 
them using the e-expansion. This certainly requires some assumptions about the underlying 
physics. 

As has been mentioned in the Introduction, a unique property of the critical states is mul- 
tifractality. This property is characterized by the power-law dependences of averaged powers 
of eigenfunction amplitudes |\l/£;(r)|. Two of such power-law dependences are knownll'ii. 
One of them determines the scaling of a single eigenfunction with respect to the size of the 
system L. 

J2{\^4r)f^{E - E^)) oc L-'^'^^-^), (4.5) 

r.n 

Another one determines the correlations of different eigenfunctions as a function of energy 
difference: 

Y: (l^n(r)r \^m{r)ME - E^) 6{E' - EJ) ^ \E - i^'^ ^'^"^^ (4-6) 

r,n,m 

In Eqs.( [4.5|) , (|4.6| ) dg < d is a. fractal dimension that depends on q {^^multifractalit'i/^). 
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It is remarkable that Eq.( ^4.5| ) can be derivedEl for the case of 2d metals by means of 
an instanton approximation similar to the one we used in this paper. The spectrum of the 
fractal dimensions dq obtained in this approximation turns out to be linearH: 

Tj 2 

dq = d--q, r] = d-d2 = —, (4.7) 

where P = 1,2,4 for the orthogonal, unitary and symplectic ensembles. 

It is reasonable to assume that the power-law tail in P{X) is another signature of 
multifractalityS. Then one may hope that the expression for n in terms of the struc- 
tural constant of multifractality t] = d — d2 providesa better approximation for /i than the 
e-expansion. By using Eq.( |4.71 ) and the relationships^ ^ = 2l3g2, we obtain: 

/i = -. (4.8) 

V 

The derivation of Eq.( [4.8|) that we have just carried out for the d = 2 case is based on two 
crucial facts: i). the exponent fi is determined by the spectrum of multifractality dq and ii). 
this spectrum is linear (for q <^ We will now make the assumption that i). is valid 

for any critical state. Since for any critical state with weak multifractality the spectrum 
of dq is expected to be linear up to very large values of q, we believe that Eq.( [4.8|) is valid 
for any critical state with weak multifractality. In contrast to Eq. ( [4.3|) , the relationship 
between n and t] Eq. ([4.8|) is independent of dimensionality and the symmetry parameter 
j3 and should apply to M critical states in the Quantum Hall regime and for systems with 
spin-orbit interaction!^^. 



V. NUMERICAL EVALUATION OF P{K) FOR THE ANDERSON MODEL. 

For our numerical analysis we consider a tight-binding model on a square lattice of L'^ 
sites. The one-particle Hamiltonian is: 

H = Y^ e4c, + t Yl ie'''' c\c, + e"^^- c]c,). (5.1) 

i {<ij)> 

The site energies are randomly distributed with uniform probability between —W/2 and 
W/2 in units of t = 1. The parameter W controls the amount of disorder in the system. 
The phase shifts 6ij in the hopping term connecting nearest neighbors represent the effect of 
an external perturbation that breaks the T-invariance of the system. As for the analytical 
calculations we consider two types of such perturbations. The first one (case I) is the usual 
Aharonov-Bohm flux $ = (0/27r)$o that pierces the system closed to a ring geometry giving 
rise to a global shift of the boundary conditions in one direction. We will choose a gauge such 
that each hop in x-direction picks up a phase 9ij = (p/L, so that total twist of the boundary 
condition is (j). The second one (case II) is a random magnetic flux. In this case the gauge 
is such that the phase 6ij relative to a hop in the x-direction is Gaussian distributed with 
zero average and variance equal to {9fj) = (0/L)^. For this gauge, the vector-potential 
A(i) oc {9ii^i,0,...} is defined on the dual lattice with sites in the middle of bonds in 
x-direction. Thus this is a random vector-potential model with a short-range correlator 
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{Aa{i)Aj3{j)) = {^/ Sa,xSi3,x Sij of the type given in the continuous approximation by 
Eq.( p.6|) . The only difference is that the correlator is anisotropic. This difference is not 
important, since it leads only to a constant factor 1/d in Vr that can be absorbed in the 
parameter 0. This kind of perturbation is qualitatively different from case I, since it acts 
locally. 

The numerical evaluation of the curvature is based on the representation of the second 
derivative by the finite difference: 

= 2<p~^EM - EniO)]- (5.2) 

In using this formula one should take care that is small enough in order for Eq. ( p.2| ) to 
be valid. On the other hand, a too small would result in big numerical errors because of 
the finite numerical precision in evaluating En. The optimal choice of should be made 
for each level En separately, since the level curvatures vary in a wide range for a given 
realization of disorder. To this end, we diagonalize the Hamiltonian for several values of 
(up to ten values) for each disorder realization and choose the smallest for which the 
normalized difference |-En(0) — En{0)\/ En{0) is still larger than some conveniently chosen 
small parameter. In order to attain a smooth curve and decrease the statistical fluctuations 
the statistical average has been made both over the energies (in the energy window of width 
4 centered ai E = 0) and over many realizations of the disorder (there were typically few 
thousands of them). 

A. Corrections to P{K) beyond RMT in 3d metals. 

In this section we compute the finite g corrections to the shape of the curvature distri- 
bution in the metallic regime, comparing the numerical results with the regular corrections, 
Eq.(^,(^). 

It turns out that the magnitude of corrections is small and we need to consider a rather 
large disorder {W ~ 10) to detect it. At yet larger values of disorder we may enter the 
critical regime. The onset of the critical regime exhibits itself in the weak dependence of the 
conductance g and the magnitude of correction to P{k) on the system size L. In contrast 
in a good metal a naive estimation of Cd in Eq.(|2.1|) yields |Crf| oc W^/L"^. The upper limit 
of the coefficient |Crf| can be estimated from the critical conductance g'^ = 4:7r^g* ^ 1. Since 
\Cd\ is proportional to a small parameter {47T^g)~'^ = {ng^)""^, its value just near the critical 
region is proportional to a small numerical factor I/tt^ ~ 0.1. It is this small numerical factor 
together with the strong ly-dependence of that makes the correction to P{k) small and 
Eq.( p.l| ) qualitatively applicable even very close to the mobility edge. 

1. Global Vector Potential. 

In Fig. ^ we show the numerical results for 6P{k) = P{k) — Pwoik) for the 3(i metallic 
regime in case I. The calculations are performed for system size L = 8 and disorder W = 12. 
The number of disorder realizations is 1500. The deviation of P{k) from the RMT result 
is very small, less than one percent. The magnitude of the statistical noise present after 
averaging is done appears to be only a little smaller than the signal itself. Nevertheless the 
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general trend of the curve agrees with the analytical prediction Eg. ( p.l| ): P{k) is above the 



RMT result at small k. We have used the coefficient Cd in Eq.(2J.) as a free parameter in 
the least square ffiting of the numerical results. The value C3 = 0.0044 found from such a 
fitting is probably a reliable estimate of the magnitude of the correction SP{k) in the above 
case. 

2. Random Magnetic Flux 

The same correction 6P{k) for the case of a random magnetic flux is displayed in Fig. |^. 
The values of the parameters of the Hamiltonian are the same as for the previous case. 
Despite the statistical fluctuations are still rather strong, the numerical results are quite 
significant. We see that again the expression Eq.( p.l|) provides a rather good one-parameter 
fitting function for the numerical results. However, in this case the coefficient = —0.014 
is negative in full agreement with the analytical prediction. Moreover, the numerical results 
are consistent with the analytical prediction even quantitatively. It follows from Eg. ( p. 2]) 
that there is a magic relationship for the ratio of amplitudes of the correction in case I and 
case II: 

R=^ = --, (5.3) 
^(11) 3' V ; 

Our calculations give a result R = —0.32 which is in an amazingly good agreement with 
Eq.(U. 

B. P{k) at the mobility edge in 3d and in 2d metals. 

A numerical investigation of the distribution P{k) at the Anderson transition critical 
point has been already carried out in Ref . [ ^ . The main finding of the numerical simulation 
is that the distribution function at the mobility edge is remarkably well fitted by the formula: 

P^(k) = — -. (5.4) 

(l + |A;|(2-"))2^ ' 



with a ~ 0.4. Equation ( |5.4|) defines a function that has a branching point of the type 
Eq.( [4.4D at = and the asymptotic behavior Pa{k) oc \k\~^, expected in all cases where 
there is energy level repulsion R{uj, 0) oc |a;| at a; -C 1. The function Pa{k) is a rather special 
one0, since once it is properly normalized by choosing Aa = (2 — a)r[3/(2 — a)]{[r[l/(2 — 
a)]r(2/(2 — a)]}""*^ > 1, it automatically satisfies the condition / dk\k\Pa{k) = 1 for all a. 
Thus the full distribution Pa{k) is determined uniquely by its value at k = 0. If we take 
Pa{k = 0) equal to the numerical result (that is known up to a small error bar) the whole 
curve is parameter free. Alternatively we can consider a as a free parameter that should be 
determined by a least square fitting of the overall numerical curve. These two procedures 
yield very close values for a and an excellent overall fitting of the numerical curve, implying 
that the extrapolation by means of the function Pa{k) is very self-consistent. 

In Fig. P we plot the results for the difference 6P{k) = P{k) — Pwo{k) for the critical 
disorder W = 16.5 and the system size L = 12 as compared to two one-parameter fitting 
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curves provided by Eg. (p.lD and Eq. (|5.4] ). It is clearly seen that despite the analytic function 
(5Preg(fc) given by Eq.( p.l|) reproduces a correct qualitative behavior of 6P{k), there is some 
feature at small \k\ that is captured better by the nonanalytic fitting function, Eq.( |5.4| ). 

We note also that this enhancement at small \k\ relative to (5Preg(^) is size- dependent, 
with P{k = 0) increasing with L so that a full saturation is not reached even for L = 10. Such 
a behavior is beyond the one-parameter scaling and can be explained by the termination 
of the power-law tail in the characteristic function P{X) at a finite A = {L/l)^. Indeed, 
the power-law tail given by Eq.( |4.3| ), makes a contribution to P{k = 0) proportional to 
j^C^-ZO dx This contribution is size-dependent and increases with increasing L (with 
a saturation at L/l —>■ oo) even exactly at the critical point where the exponent /i is a 
constant. 

Now let us calculate the fractal dimension d2 = 3 — t] using the relationship, Eq. (|4.8| ), 
between rj and n = 3 — a. For a ~ 0.4 we have fi ~ 2.6 and d2 ~ 1.5. This value is in a 
good agreement with direct evaluation^ of d2 from Eqs.(^l5| ) ,( |4.6| ). 

We also checked that in a 2d metal the deviation from the RMT result 6P{k) has the 
same qualitative form as in the 3d critical case. In Fig. |^ we present plots of 6P{k) for 
different values of L and W = 6. An interesting feature seen in the figure is the fixed point 
in 6P{k) a.t k ^ 0.35. Despite the one-parameter function Pa{k) does not have an exact 
fixed point at the same value of fc, all the curves obtained varying a do get very close to 
each other at /c ~ 0.35 in what looks almost like a fixed point. As in the 3d critical case, 
the analytic function 6P^cg{k) fits well the overall distribution but fails to describe the sharp 
enhancement for small curvatures that is instead well described by the function Pa{k). This 
is shown in Fig. ^ for the system L = 30, = 6. 

VI. CONCLUSIONS 

In this paper we have investigated both analytically and numerically the relationship 
between the statistics of eigenf unctions and the spectral statistics in disordered conductors. 
The level curvature distribution has been chosen as the target of our study, since it is the 
simplest known example of parametric spectral statistics that can be used as a spectral 
probe of the structure of eigenfunctions. 

The main results of the paper are formulated in Eqs.(|2.1|),(|2.2|), and Eqs.(|4.4]),([4.8D. Numer- 
ical results in agreement with these analytical predictions are presented in Figs. ^J^, and 
respectively. The first two equations describe the regular corrections beyond RMT to the 
level curvature distribution in disordered metals. These corrections stem from long-range 
correlations in the wave functions with a typical length scale of the order of the sample size. 
The latter two equations summarize the effect of local irregularities (sharp peaks) in the 
structure of eigenfunctions in its most developed form (multifractality) in the critical region 
near the mobility edge. 

Moreover, Eq.( [4.8|) suggests an explicit relationship between the fractal dimension 
d2 = d — f] of a critical eigenfunction and the exponent /i in the power-law tail Eq. ( |4.3D 
of the characteristic function -Pc(A) that describes the particular parametric spectral statis- 
tics considered here (level curvature distribution). 

Note that Eq. ( |4.^ ) is more general than Eq. ( |4.4|) . The latter requires rather strong 
multifractality > |, while the former applies to a generic critical state. For instance, it 
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would be interesting to check its validity for the critical state in the quantum Hall effect, 
where rj 0.5 and we predict /i ~ 8. Recent progressill in numerical simulations on the 
Chalker-Coddington network mo deli seems to make the task attainable. 

As far as the regular corrections are concerned, there is an interesting question of what 
happens to them for d > A. The sum in Eq.( p.2| ) is a parameter-free number only for 
d < 4 when it converges. For d > 4 the sum is divergent and requires a cut-off at large 
|q|. Thus for the correct evaluation of this sum it is necessary to go beyond the diffusion 
approximation and the approximation of slow spatial variations of the field Q{r) in the 
nonlinear sigma- model. The divergent sum in Eq.( |2.2| ) implies that for d > 4 the correction 
6P{k) to the level curvature distribution is dominated by short-range spatial correlations of 
the eigenfunctions, in contrast to d < 4 where it is dominated by long-range correlations. 
Based on our discussion at the end of Sec. |11 A|, the short-range nature of the eigenfunction 



correlations in c? > 4 is most likely the cause of the change of sign in the correction 6P{k) 
for the "global" case I (as compared to d < 4) and the reason why 5P{k) shows qualitatively 
the same behavior of case II ("local" perturbations). 

One of the most important results Eq.( |2.4| ) of our calculations is that the ratio r{g) of 
the mean level curvature {\K\) and the mean Drude conductance 2g is not a constant and 
is always larger than the RMT result r = 1. This is in a qualitative agreement with the 
resultaS that {\K\) oc ^ in a strictly one-dimensional case where (? -C 1 and localization 
effects are strong. Indeed, let us assume that the square-root dependence is typical for 
strongly localized states in any dimensions. Then the function r{g) should behave like 
r{g) oc g'^/"^ at small g and r{g) — > 1 for g — > oo. If in addition we make the natural 
assumption that r{g) is a monotonic function, we arrive at the conclusion that 6r{g) > 



everywhere in agreement with Eq.(2.4). 

A similar deviation from the proportionality relationship {\K\) oc g has been observed 
recently in numerical simulationscj. 
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APPENDIX A: EFFECTIVE ACTION 



The parametrization ( p.l3|) enables us to single out fast modes Q in the action ( p.l2| ) as 
follows: 

f[q,Qo] =5tr|^(P^,g] - tLWQf - (t^k<l)'QA + Qj^ qY (A1) 

where the following notation has been used: 

= To$T,-\ Qa = ToAT,-\ Qj = T^JT^K 
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g = ^ = uDL''-\ STr...^ J yStr... , K = 2g ■ k. 

Taking advantage of the fact that we are interested in hmit only, it is convenient to 
absorb the large conductance into the definition of the flux ^g<f)'^ — > 0^ (wherever possible). 



Then introducing the dimensionless gradient d = i/7r/2LV and using the identity 



-STrQAQ = STrQ^Q, 



we recast the action in the following form 
F 



Q,Qo\ =STr{j{['g^,Q]-z^dQy -Uk^l + Qj]Q\. 



(A2) 



These manipulations, along with fact that W scales asW = g ^^"^w, w oc 1 enable us to con- 
struct a perturbative expansion in the parameter 1/(7— straightforwardly, simply expanding 
F Q, Qq in powers of w. This procedure is equivalent to selecting diagrams contributing 
to the same order in 1/g. Then for the partition function we get the representation: 

a2 



lim 



*o 327Tg 



3? / DQoDQ ■ J Q 



exp 



(A3) 



where J 
Appendix 



Q 



is the Jacobian of the transformation Q —>■ (^Q, Qo) (obtained in the following 
). The action (K3) is expanded as follows: 



-| 00 
Q, Qo] = F [Qo] + -STr {dwf + Q-'^'^F^,, [Qo, 

^ n=l 



(A4) 



The first term of the expansion is nothing but the zero- mode action (i?MT-limit): 



F[Q, 



STr 



1 /— > \2 



The second term in ( |A4| ) corresponds to the noninteracting diffusion modes approximation. 
The other terms in ( [A^ ) describe the interaction of the diffusion and zero modes with g~^ 
playing the role of a coupling constant. For our purposes it is enough to keep the first four 
interaction terms of expansion: 



Fi/2 = iSTrQ ^wdw, 



Fi = ^STr i^AAw^ + ^Awf + ^ (dwfw^^ + 
+ ^.STr {"^^ {dw^ -2wdww)}, 

F3/2 = \sTr |aAm;3 ^ 2 (^'^ ^Awf + 



(A5) 



(A6) 



(AT) 
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(A8) 



We have used the notation: 



A= Q^AQ^-ikQ^-Qj. 



(A9) 



Averaging over the fast modes can be implemented by the use of contraction rules derived 
in@. We recast them in the following form: 



'W (r) RW (r')) = H (r, r') [STrR - STrARA -R + ARA} 
'W (r) STrRW (r')) = U (r, r') [r -R - A (^R - r) A} . 

The propagator 11 (r, r') satisfies the diffusion equation: 



(AlO) 



- An fr, r') 



1 



irgL^ 



Sir 



n-'icO = ng (qL)^ 



(All) 



The contraction rules ( |A1CI| ) provide a basis for integrating out the fast modes in ( [A 31 ). 
For example, straightforward but lengthy calculations give us the following rules for the 
integration of product of two vertices from ( |A^ ) containing no gradients: 

(^STrAAW"^ (r) STrAAW^ (r')) = 4^ (r, r') {[STrAAf - [STrA]'^} , 

/ 



STr 



BAW r 



STr 



+ 



STrB' 



STrB^A 



BAW (r'; 

2 



(r, r') ■ STr (BA) 



+ 



+ STr UbA)^ + B^-2 (e^a)' 



STr 



BAW (r) 



STrAAW"^ (r" 



(r, r') ■ STr A (2B^A - BAB - A (BA)' 



where the matrices satisfy the symmetry relations A = A, B = —B (for notations seeQ)and 
the brackets (,,,)^ mean keeping only the connected parts of the correlators after averaging 
over fast modes W . 

These calculations can be easily extended to couplings of gradient vertices as well. Then 
having at hand all possible couplings we are able to construct a cumulant expansion for the 
partition function ( |A3| ) . Keeping all nonvanishing correlators up to the second order in 1/g 
we arrive at: 



F^" = Fo + {F,)-l{F^)^ + l{F, 
1 



Fl/2 



4! 



1/2 



(A12) 



eff 



(Fi) 



J2U{cO-STr(Q^A 



2 (^1^1/2 



4^n2 (q).(/2 + /3 + /4). 
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Here the term /2 corresponds to coupling of two gradientless vertices ( |A7|) : 



-8/2 = [STrAXf - [STrAf + 



STrQ.A 



STr{{Q,A) -2{Q,A] + 2A {2Q ,A - Q ^AQ ^ - A {Q ,A 



/s describes the couphng of two gradient vertices ( |A5| ) to the gradientless one(^) 



2d- h = STrAA ■ STr (^^A)^ - 2STrA ■ STrl^^A + STr (^^A 
-4 



n 2 



STt'^Ia +2STr('^lA 



and the last contribution /4 comes from 4 coupled gradient vertices 



9 



d{d + 2)- U 



STr(Q.A 



1 2 



4 



STrQ.A 



The expression for F^-^-^ can be simplified considerably if one uses the following factor- 
ization properties of long traces: 



Str[{4>Qo)'] 



(A13) 



5tr[(0go)'02Qo] = ^5tr[(0go)1 Str[cf'Q,] = ^ 



a 

Y 

2i oj-»,rl2/ 



2 



a 



The factorization holds in the leading power in the noncompact angles in the original Efetov's 
parametrization and can be shown by straightforward but extremely lengthy calculations. 

It is worth noting that the Jacobian J Q also contributes to the effective action (the 
last term in Eg. ( p.l9|) ). Taking into account Eq. (|B10D of Appendix |^ one concludes that 
the Jacobian leads to the replacement + ^ in n(q). Thus the l/g"^ contribution from the 

Jacobian follows from the correction to 11 (q) in F^-^^: 



(A14) 



Collecting all the results obtained in this Appendix we arrive at the effective action JF [Qq] 
given in Eqs.(^,(CT)- 

Using these equations one can represent P{k) in the differential form of Eq.( |2.22D . It 
turns out that all the relevant derivatives of Pa{k) over k can be expressed through the 
derivatives over a. Using the following identities: 
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1 + 



a=l 



d_ 

da 



Pa{k) 



d 



(A15) 



a=l 



+ Pa{k) 



a=l 



da? da 



a=l 



d^ 



Pa{k) 



d d^ 



a=l 



da da^ 



Pa{k) 



a=l 



k 



dk'' 



■Pa{k) 



d d^ 

2 + 4— + 



a=l 



da da^ 



Pa{k) 



a=l 



we can finally represent 5P{k) in the form Eq. (|2.24| ) that contains only the first and second 
derivatives with respect of a. 



APPENDIX B: CALCULATION OF THE JACOBIAN 

In the evaluation of the Jacobian of the transformation Q — > (^Qo, we follow the 
procedure proposed in Ref. [ ^ and prove that the Jacobian does not depend on zero-mode 
Q(^. The derivation is simplified if we go to the "rational" parametrization first: 

Q = {l~W)A{l-W)-\ W=(^^^y (Bl) 

This parametrization has been known since Efetov's worki. The advantage of this represen- 
tation lies in the fact that Jacobian of the transformation Q ^ W equals to unity. On the 
another hand we have the decomposition ( |2.13| ). Taking into account that Tq belongs to the 
graded coset space UOSP (2, 2|4) /UOSP (2|2) ® UOSP (2|2) the complete parametrization 
takes the form: 



Q = T,''QTo, To = ^^^, Q={l-W)A{l-Wy\ (B2) 
Bo\ ( B 



Comparing these two representations ( pi|) and ( p2|) we derive the connection between W 
and (Wq,W): 



W = {Wo + w){1 + WovjY\ w={i-W'^) ^'^W{l-W^)^'^, (B3) 
B = (fio + h) (1 + B^h)-^ , b=(l- B^bX^'^ b(1- BoBo, 



Since the field b (being proportional to W) contains only nonzero momenta it can be treated 
perturbatively. Going to the Fourier representation and expanding ( [B3| ) up to second order 
in b we separate zero- and nonzero momenta: 

5 (k = 0) = 5o - 56qfio&-q, (B4) 

B(k^O) = S (b^- 6k+q5o^-q + ^k+qi+qa^O^-qi^O^-qa) 5 
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where S = (l - 5oBo) has been introduced and the summation over repeated indecis is 
imphed. Therefore, we are interested in Jacobian of the transformation B — > {Bq, b) which 
is equivalent to (-B (k = 0) , -B (k 7^ 0)) {Bq, Bi^^q). The corresponding Jacobian can be 
represented in the block-matrix structure: 



J = SDet 



dB (0) /dBo dB (0) /dh 



(B5) 



dB (0) /BBq dB (0) /dh j ' 

where a short notation for the following supermatrices has been introduced: 

[dB{Q)/dh],, =dBiq = 0)/dh, 
[dB{(l})/dBoU =dB{kj^O)/dBo, 
[a5(0)/96],,, =dBik^O)/dbk', 

and right derivatives are impliedi. Using the identity for the superdeterminant of block 
matrices we recast ( P5|) in the following form 



J = Jl-J2 



Jo = S'det 



Ji = SDet 



dB (0) 



(B6) 



dB (0) dB (0) fdB (0) \ dB [k' ^ 0) 



dB, 



dh 



dh 



kk' 



dB, 



where S'det acts within the space of 4 x 4 matrices, while SDet spans fc-space also. Lengthy 
but straightforward calculations give us following very useful formulae: 



d (DbF) 
Str \^ ' = StrD ■ StrF, 
do 



Str 



d (DhP) 



dh 



StrDF, 



Str 



d (DihFi) d {D2b2F2) 



dbi 



dbo 



StrDiD2 ■ StrFiF2, 



d [DihiFi] d [D2h2F2 ^ 

Str ^ I ' ^ = StrDiD2F,F2, 



dbi 



dbo 



(B7) 



Str 



dbi 



d {Gb) 
db 



^d(D2b2F2 



dbo 



Now we expand Ji up to second order in b: 



where 



K = exp < STr 



diSb)\-'d{bBobBob) 1 



StrDiG-^D2FiF2. 



db 



J, . soJ-m . K. 

db 



STr 



db 



d{Sb)\-'d[SbBob 



db 



db 
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Taking into account that SDet \d (Sb) /db] = 1 and exploiting 

'd{Sb)\-'d[SbBobBob 



we obtain: 



Str 



Str 



db 

d{Sb) \ 
db J 



db 



d SbBnb 



db 



J2 StrbqBoStrBob^^, 

q^O,k 



2 Strb]^>^]^BoStrBobk_k' 

kYO.kT^O 



Finally, we arrive at: 



Ji = exp 1^ 5'tr6_ki?o5'tri?o&k| • 



Then with the same accuracy we may recast J2 in the following form: 
J2 = 5det(l -5J2), 

5J2 



diSb^Bob^^] diSb^Bob^^] fd{Sb)\-'d{Sb^) 



In J9 Str 



dBo 



db. 



db 



dBa 



dBo 

Using formulae ( [B7[ ) we obtain: 

d [Sb^Bob^ 



-Str- 



d[sb^Bob^^) f d iSb)Y' d {Sbi, 



db. 



db 



dBo 



Str- 



Str 



dBo 

diSb^Bob^^ 



StrSbqb^q — StrBobq^Bob^f^ 



'd{Sb)Y^ d{Sb 



db]^ \ db 
1 . 



dBo 



= -StrSb^k [1 - BoBo) b^ - StrBob^StrBob^k 
+ StrSb^]J)i^ — StrBob\f_Bob^]^ 
As a result we find for the second contribution J2: 

- In J2 ^ Str (1 - BoBo) 6-k (l - ^oSq) K + StrBob^StrBob-^ 



dr 
V 



StrW^ + StrBobkStrBob^k 



(B8) 



(B9) 



The parametrization used in this Appendix differs slightly from the one used in the main 
body of paper. To come back to the original parametrization we must substitute W ^ W /2 
. Then collecting contributions from Ji and J2 together we finally obtain Jacobian in the 
form: 



J 



Q 



exp 



dr 

8J V 



StrW 



(BIO) 



As one can see from Eq. (piOp the Jacobian does not contain the zero-modes at all. Morover, 
being quadratic, it just plays the role of a small frequency (~ 1/g) in the free diffusion 
propagator. Expanding such modified propagator we recover the Jacobian contribution to 
the effective action Eq. (|A14|) . 
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APPENDIX C: SOLUTION TO 2D LIOUVILLE EQUATION 



On the boundary of the square Vt [see Fig. |^ the function u{z) and its first derivative 
must be continuous. Because, by contruction, u{x) is anti-symmetric when reflected around 
the sides of fi, it must be zero on the boundary of Vl: 



dF 



dz 



= ^{i+\F{z)\r, (ci) 



For the case 7 — > we look for a solution ^ 1. Then we have 

d f I \ 7 



dz \F{z) J 4 



exp(ze), ImQ = 0. (C2) 



Suppose we manage to find Q{z) which is an analytic function inside the square Q and is 
real on its boundary. Then the equation for 1/F is trivially reducible to quadratures and 
the solution F{z) is an analytic function inside the square Q. Then in the limit 7 <^ 1 the 
solution to the Liouville equation that obeys the condition u{z) = on the boundary of the 
square fl is found by substituting the function F{z) into Eq.( p.26| ). In particular, in the 
region where ^ 1 we have: 



Consider the function 



which does the conformal transformation of the unit circle in the complex plane of t onto 
the square fl in the complex plane of z. This function obeys the symmetry property: 

z{it)=iz{t), z{t*)=z*{t). (C5) 

If we choose e*®*^^-* = where A; is a real parameter, the analytic function Q(z) will 

automatically be real on the boundary of the square fl. Integrating Eq. (|C^ ) we arrive at: 

1 



Fiz) 

For |z| <^ 1 we have from Eq.( 



and 



7 Vv^ 

We see that, indeed, ^ 1 for 



lJ\t{z'r'dz'. (C6) 



— t (C7) 
vr 



Fiz) = --(^] {k-l)z''-\ (C8) 
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7 



(A;>1). 



vr [2nV2{k-l)_ 
In this region the function u{z) is given by: 

u{z) = -2kln\t{z)\ = -2kmnt{z), {\z\ > tq). 



(C9) 



(CIO) 



Thus for ro <^ 1, the solution to the 2d Liouville equation with boundary conditions 
u{z) = on the boundary of the square fl is given by Eq.( p.26|) with F{z) defined in Eq . (|C^ ) . 
For 1^1 <^ 1 this solution for u{z) depends only on \z\ = r and is given by Eq.( 
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FIGURES 



FIG. 1. Perturbative 1/g correction to the curvature distribution P{k) beyond the RMT result 
-Pwd(^)- The expression for (5Preg(fe) = P{k) — PwD{k) is given in Eqs. ( |2.1D ,(|2^). The dashed 
hne represents the case of global T-breaking perturbation (Case I). The solid line represents the 
case of local T-breaking perturbation (Case II). 

FIG. 2. Plot of the periodic solution u{x) of the Liouville equation Eq. ( |3.23| ), that enters 
the instanton action of Eq. ( |3.22| ), for the quasi-lD case. The function is defined on the ring 
— ^ < X < ^. In the interval —j < x < j the function u{x) is the positive solution of Eq. ( ^.27 ) 



with the choice 6=0 and k found from the condition n(ibl/4) = 0. In the intervals ib^ < x < ib| 
u{x) is constructed by antisymmetric reflection around the points ±4. 



FIG. 3. Definition domain of the periodic solution u{z) of the Liouville equation, Eq. ( p. 23 ) 



for the 2d case. The domain is a torus —1/2 < x,y < 1/2, is represented here by the large square. 
Inside the square with vertices at z = ±1/2, ibi/2 the function u{z) is the real positive solution 
of Eq. (|3l2t ). In the remaining part of the larger square, u{z) is constructed by anti-symmetric 
reflection around the sides of the square and its sign is negative. 

FIG. 4. Numerical results for 6P{k) = P(A;) — Pvvd(^) for the 3d Anderson model in the metallic 
regime. The disorder isw = 12t and the system size is L=12. The curvatures are calculated for case 
I, global vector potential. The smooth curve is fit with the analytical result given in Eqs. (p.l|),( ^l2| ). 
The coefficient C3 is taken as a free parameter in the least square fitting. 



FIG. 5. The same as in Fig. |^ but for case II, local random magnetic field. 

FIG. 6. Numerical results for the difference 6P{k) = P{k) — PwD{k) for the 3d Anderson model 
of system size L = 12 at the Anderson critical point (disorder W = 16.5). The dashed and solid 



curves are two one-parameter fitting curves provided by by Eq.(2.1) and Eq.( 



FIG. 7. Numerical results for the difference 5P{k) = P{k) — P\Y£){k) for the 2d Anderson model 
of different system sizes L = 12, 16, 20, 24, 30 and site disorder W = 6. 

FIG. 8. Numerical results for the difference 6P{k) = P{k) — Pwoik) for the 2d Anderson model 
of system size L = 30 and site disorder W = 6. The dashed and solid curves are two one-parameter 



fitting curves provided by by Eq.(^.lD and Eq.(5.4). 
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